#include "MUltiPhasePyrolysisMaterial.h"
#include "ComputeRhoMaterial.h"
#include<iostream>
#include "libmesh/quadrature.h"
using std::cout;
using std::endl;

registerMooseObject("TrilobitaApp", MUltiPhasePyrolysisMaterial);

template<>
InputParameters validParams<MUltiPhasePyrolysisMaterial>()
{
  InputParameters params = validParams<Material>();
  //add regular thermal properties
  params.addParam<Real>("surface_emissivity", 0, "emissivity");
  params.addParam<Real>("R", 287, "gas constant");
  params.addRequiredParam<std::vector<Real> >("kv_list", "The vector of kv values of origin material for building the piecewise function");
  params.addRequiredParam<std::vector<Real> >("kvT_list", "The vector of temperature values for building the piecewise function to compute kv ");
  params.addRequiredParam<std::vector<Real> >("kc_list", "The vector of kc values for building the piecewise function to compute kc ");
  params.addRequiredParam<std::vector<Real> >("kcT_list", "The vector of temperature values for building the piecewise function to compute kc ");
  params.addRequiredParam<std::vector<Real> >("cpv_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
  params.addRequiredParam<std::vector<Real> >("cpvT_list", "The vector of cpv values for building the piecewise function to compute thermal conductivity of origin material ");
  params.addRequiredParam<std::vector<Real> >("cpc_list", "The vector of cpc values for building the piecewise function to compute cpc ");
  params.addRequiredParam<std::vector<Real> >("cpcT_list", "The vector of temperature values for building the piecewise function to compute cpc");
  params.addRequiredParam<std::vector<Real> >("hv_list", "The vector of cpv values for building the piecewise function to compute hv ");
  params.addRequiredParam<std::vector<Real> >("hvT_list", "The vector of cpv values for building the piecewise function to compute hv ");
  params.addRequiredParam<std::vector<Real> >("hc_list", "The vector of cpc values for building the piecewise function to compute hc ");
  params.addRequiredParam<std::vector<Real> >("hcT_list", "The vector of temperature values for building the piecewise function to compute hc");
  params.addRequiredParam<Real>("rhov","density of original material");
  params.addRequiredParam<Real>("rhoc","density of fully charred material");
  //add special parameters to compute pyrolysis Arrhenius relation
  params.addRequiredParam<Real>("pyrolysis_heat","pyrolysis heat");
  //add coupled variables
  params.addRequiredCoupledVar("temperature", "Coupled Temperature");
  //add porous media parameters
  params.addRequiredParam<Real>("porosity_v", "porosity of original material");
  params.addRequiredParam<Real>("porosity_c", "porosity of fully charred material");
  params.addRequiredParam<Real>("permeability_v", "permeability of original material");
  params.addRequiredParam<Real>("permeability_c", "permeability of fully charred material");

  return params;
}

MUltiPhasePyrolysisMaterial::MUltiPhasePyrolysisMaterial(const InputParameters & parameters) :
    	    Material(parameters),
//		_property_pack(declareProperty<PropertyPack>("propertypack")),
		_epsilon(declareProperty<Real>("surface_emissivity")),
		_sigma(declareProperty<Real>("sigma")),
		_R(declareProperty<Real>("gas_constant")),
		_kv(declareProperty<Real>("kv")),
		_kc(declareProperty<Real>("kc")),
		_ks(declareProperty<Real>("ks")),
		_cpv(declareProperty<Real>("cpv")),
		_cpc(declareProperty<Real>("cpc")),
		_cps(declareProperty<Real>("cps")),
		_hv(declareProperty<Real>("hv")),
		_hc(declareProperty<Real>("hc")),
		_hs(declareProperty<Real>("hs")),
		_rhov(declareProperty<Real>("rhov")),
		_rhoc(declareProperty<Real>("rhoc")),
		_pyrolysis_heat(declareProperty<Real>("pyrolysis_heat")),
		_porosity_v(declareProperty<Real>("porosity_v")),
		_porosity_c(declareProperty<Real>("porosity_c")),
		_porosity(declareProperty<Real>("porosity")),
		_permeability_v(declareProperty<Real>("permeability_v")),
		_permeability_c(declareProperty<Real>("permeability_c")),
		_permeability(declareProperty<Real>("permeability")),
		_char_percent(declareProperty<Real>("char_percent")),
		_rhos(getMaterialProperty<Real>("solid_density")),
		_T_value(coupledValue("temperature")),
		_gradient_T_value(coupledGradient("temperature"))
{
	_epsilon_value = getParam<Real>("surface_emissivity");
	_R_value = getParam<Real>("R");
	_kvT_list = getParam<std::vector<Real>>("kvT_list");
	_kv_list = getParam<std::vector<Real>>("kv_list");
	_kcT_list = getParam<std::vector<Real>>("kcT_list");
	_kc_list = getParam<std::vector<Real>>("kc_list");
	_cpvT_list = getParam<std::vector<Real>>("cpvT_list");
	_cpv_list = getParam<std::vector<Real>>("cpv_list");
	_cpcT_list = getParam<std::vector<Real>>("cpcT_list");
	_cpc_list = getParam<std::vector<Real>>("cpc_list");
	_hvT_list = getParam<std::vector<Real>>("hvT_list");
	_hv_list = getParam<std::vector<Real>>("hv_list");
	_hcT_list = getParam<std::vector<Real>>("hcT_list");
	_hc_list = getParam<std::vector<Real>>("hc_list");
	_func_kv_T = LinearInterpolation(_kvT_list,_kv_list);
	_func_kc_T = LinearInterpolation(_kcT_list,_kc_list);
	_func_cpv_T = LinearInterpolation(_cpvT_list,_cpv_list);
	_func_cpc_T = LinearInterpolation(_cpcT_list,_cpc_list);
	_func_hv_T = LinearInterpolation(_hvT_list,_hv_list);
	_func_hc_T = LinearInterpolation(_hcT_list,_hc_list);
	_rhov_value =  getParam<Real>("rhov");
	_rhoc_value =  getParam<Real>("rhoc");
	_pyrolysis_heat_value =  getParam<Real>("pyrolysis_heat");
	_porosity_v_value =  getParam<Real>("porosity_v");
	_porosity_c_value =  getParam<Real>("porosity_c");
	_permeability_v_value =  getParam<Real>("permeability_v");
	_permeability_c_value =  getParam<Real>("permeability_c");
}
void MUltiPhasePyrolysisMaterial::initQpStatefulProperties()
{
	//compute regular thermal properties for orinin and charred material separately
	_epsilon[_qp] = _epsilon_value;
	_R[_qp] = _R_value;
	_kv[_qp] = _func_kv_T.sample(_T_value[_qp]);
	_kc[_qp] = _func_kc_T.sample(_T_value[_qp]);
	_cpv[_qp] = _func_cpv_T.sample(_T_value[_qp]);
	_cpc[_qp] = _func_cpc_T.sample(_T_value[_qp]);
	_hv[_qp] = _func_hv_T.sample(_T_value[_qp]);
	_hc[_qp] = _func_hc_T.sample(_T_value[_qp]);
	_rhov[_qp] = _rhov_value;
	_rhoc[_qp] = _rhoc_value;
	//compute porous media parameters for orinin and charred material separately
	_porosity_v[_qp] = _porosity_v_value;
	_porosity_c[_qp] = _porosity_c_value;
	_permeability_v[_qp] = _permeability_v_value;
	_permeability_c[_qp] = _permeability_c_value;
	//compute char percent :1-rhov/(rhov-rhoc+1e-8)*(1-rhoc/rho)
	_char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
	//compute regular thermal properties for solid material
	_ks[_qp] = _kv[_qp] * (1 - _char_percent[_qp]) + _kc[_qp] * _char_percent[_qp];
	_cps[_qp] = _cpv[_qp] * (1 - _char_percent[_qp]) + _cpc[_qp] * _char_percent[_qp];
	_hs[_qp] = (_hv[_qp] * _rhov[_qp] - _hc[_qp] * _rhoc[_qp]) / (_rhov[_qp] - _rhoc[_qp] + 1e-8);
	//compute porous media parameters for solid material
	_porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
	_permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
	_pyrolysis_heat[_qp] = _pyrolysis_heat_value;
}



void MUltiPhasePyrolysisMaterial::computeQpProperties()
{
	//compute regular thermal properties for orinin and charred material separately
	_epsilon[_qp] = _epsilon_value;
	_R[_qp] = _R_value;
	_kv[_qp] = _func_kv_T.sample(_T_value[_qp]);
	_kc[_qp] = _func_kc_T.sample(_T_value[_qp]);
	_cpv[_qp] = _func_cpv_T.sample(_T_value[_qp]);
	_cpc[_qp] = _func_cpc_T.sample(_T_value[_qp]);
	_hv[_qp] = _func_hv_T.sample(_T_value[_qp]);
	_hc[_qp] = _func_hc_T.sample(_T_value[_qp]);
	_rhov[_qp] = _rhov_value;
	_rhoc[_qp] = _rhoc_value;
	//compute porous media parameters for orinin and charred material separately
	_porosity_v[_qp] = _porosity_v_value;
	_porosity_c[_qp] = _porosity_c_value;
	_permeability_v[_qp] = _permeability_v_value;
	_permeability_c[_qp] = _permeability_c_value;
	//compute char percent :1-rhov/(rhov-rhoc+1e-8)*(1-rhoc/rho)
	_char_percent[_qp] = (_rhov_value - _rhos[_qp]) / (_rhov_value - _rhoc_value + 1e-8);
	//compute regular thermal properties for solid material
	_ks[_qp] = _kv[_qp] * (1 - _char_percent[_qp]) + _kc[_qp] * _char_percent[_qp];
	_cps[_qp] = _cpv[_qp] * (1 - _char_percent[_qp]) + _cpc[_qp] * _char_percent[_qp];
	_hs[_qp] = (_hv[_qp] * _rhov[_qp] - _hc[_qp] * _rhoc[_qp]) / (_rhov[_qp] - _rhoc[_qp] + 1e-8);
	//compute porous media parameters for solid material
	_porosity[_qp] = _porosity_v[_qp] * (1 - _char_percent[_qp]) + _porosity_c[_qp] * _char_percent[_qp];
	_permeability[_qp] = _permeability_v[_qp] * (1 - _char_percent[_qp]) + _permeability_c[_qp] * _char_percent[_qp];
	_pyrolysis_heat[_qp] = _pyrolysis_heat_value;
}


